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Steady-State and Transient Boundary Element Methods 
for Coupled Heat Conduction 


Dean A. Kontinos 
Ames Research Center 


Summary 


Boundary element algorithms for the solution of steady-state and transient heat conduction are 
presented. The algorithms are designed for efficient coupling with computational fluid dynamic 
discretizations and feature piecewise linear elements with offset nodal points. The steady-state 
algori thm employs the fundamental solution approach; the integration kernels are computed analyti- 
cally based on linear shape functions, linear elements, and variably offset nodal points. The analytic 
expressions for both singular and nonsingular integrands are presented. The transient algorithm 
employs the transient fundamental solution; the temporal integration is performed analytically and 
the nonsingular spatial integration is performed numerically using Gaussian quadrature. A series 
solution to the integration is derived for the instance of a singular integrand. The boundary-only 
character of the algorithm is maintained by integrating the influence coefficients from initial time. 
Numerical results are compared to analytical solutions to verify the current boundary element 
algorithms. The steady-state and transient algorithms are numerically shown to be second-order 
accurate in space and time, respectively. 


1 Introduction 


Computational science is advancing dramatically in these first decades of the Information Age. 
Performance gains in computer technology are propelling scientific computing to the forefront of the 
engineering process. Increasing processor speed is promoting greater fidelity of the physical models, 
and increases in memory are permitting finer resolution of the physical domain. Simulations, once 
impossible to perform in a timely fashion, are routinely computed using desktop workstations. 
Furthermore, computational modeling is extending into every theatre of engineering, for example, 
fluid dynamics, acoustics, heat transfer, chemistry, astrophysics, and structural mechanics. 

Although the maturing numerical methodology is applicable across the engineering spectrum, 
much of the development is compartmentalized into the separate disciplines. For example, in the 
aerospace industry, computational fluid dynamics (CFD) simulations performed in the aerodynamics 
or propulsion group are passed, sometimes blindly, as a load condition to the structures group 
wherein a structural dynamic simulation is computed. To bridge this gap, the latest effort in compu- 
tational science is coupling analysis codes across disciplines. This effort is a natural path of develop- 
ment, but, more importantly, it is being driven by engineering systems whose complexity dictates a 
coupled analysis. For instance, in the silicon chip fabrication industry, accurate simulation of 



chemical vapor deposition and chemical etching requires modeling of the fluid dynamics in the 
reactor coupled to a surface chemistry model. In the aerospace field, advanced hypersonic concepts 
blend the propulsion system with the body mold line, thereby blurring the traditional industry 
demarcations of aerodynamics, propulsion, and structure. Also, metallic thermal protection panels 
are being implemented for the new reusable launch vehicle. Because of aerodynamic heating, the 
panels expand from the structure and alter the hypersonic flow field. Analysis of the thermal protec- 
tion system requires coupling of the aerothermodynamics to the structural response. It is in this arena 
of coupled CFD to structural analysis that the discussion is focused. 

Two basic approaches are used to solve coupled fluid/structural systems. The first is a direct 
coupling where the solution of the different fields is solved simultaneously in one large system of 
equations. Direct coupling is mostly applicable for problems where time accuracy is critical, such as 
in aeroelasticity where the time scale of the fluid motion is on the same order as the structural modal 
frequency. The second approach is a loose coupling strategy where each set of field equations is 
solved to produce boundary conditions for the other. The equations are solved in turn until an 
iterated convergence criterion is met at the fluid/solid interface. It is not within the scope of this 
paper, nor in the expertise of the author, to present a comprehensive literary review of coupled CFD- 
structural dynamic methods. Instead, reference 1 is recommended, wherein a brief bibliography and 
a coupled simulation is given. 

The loose coupling strategy is particularly attractive when coupling solid mechanics to 
hypersonic CFD. The time scales of the hypersonic flow field are often disparate from the time 
scales of the structure, thereby obviating a directly coupled time-accurate analysis. Furthermore, 
high-speed flows are replete with complex physical phenomena such as shock waves, shock-wave/ 
boundary-layer interactions, chemical reactions, and internal energy exchange. The numerics for 
hypersonic algorithms are sophisticated and, at times, temperamental. Adding structural equations 
to the system may diminish the robustness of the crafted hypersonic CFD code. Thus, the loose 
coupling strategy effectively shields the CFD code from performance degradation while increasing 
the fidelity of the global simulation. 

In the aerospace industry, finite element and finite difference methods are routinely used for the 
solution of solid heat conduction and elasticity. Well proven, these methods are readily available for 
loose coupling with a CFD code. A third option, rarely used in the aerospace community, is the 
boundary element method (BEM). Li and Kassab (refs. 2 and 3) have efficiently solved joint fluid 
and structural heating by coupling CFD to the BEM solution of the conduction in the body. The 
advantage of the BEM over a finite difference or finite element formulation is that only the boundary 
is discretized, and thus the dimensionality of the problem is reduced. As a result, it is naturally 
coupled with CFD. The boundary element grid is simply the CFD grid at the fluid/surface interface 
plus additional grid points defining the boundaries of the body. The interior of the domain is not 
discretized. This reduction in dimensionality is especially advantageous when coupling to CFD 
because interior values are superfluous; only surface values are required for coupling. Consequently, 
the BEM is potentially more efficient than finite difference or finite element methods which require 
an interior discretization to produce surface conditions. Using the BEM, the temperature can be 
computed at any desired interior point through contour integrals over the boundary solution. In 
summary, advantages of the BEM are a reduction of dimensionality, ease of discretization, and 
efficient coupling with CFD. 
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Nomenclature 


Acronyms 

BEM 

CFD 

Symbols 

A, B, C, Q 
A\, A 2 

a, b 

C p 

c 


q, c 2 

D, 


J jk ' H/j- 


// 


t55 




K 


jk 


n, n, tij 

P 

q 


4j 

Q 
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boundary element method 
computational fluid dynamics 

distance function coefficients 

integration coefficients 

nodal offset values 

domain view factor at source point 

heat capacity coefficient 

coefficients 

vector matrix containing domain integration coefficients 
internal energy per unit volume 
total energy per unit volume 
arbitrary function 

matrices containing influence coefficients 

Heaviside function 

steady-state component integrals 

transient component integrals 

tensor thermal conductivity 

scalar thermal conductivity 

length of linear element 

outward unit normal 

position vector to the point source 

vector matrix containing nodal temperature gradients 

temperature gradient boundary condition 

heat flux in index notation 

position vector to arbitrary point in the domain 

distance from source point to point in domain 

vector from source point to point in domain 
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T 

T 

A 

T 

f 

T, 

r 

T ; 

To 

t 

l f 

l 0 

u 

Uj, u k 
V 
W 

W, w 
w ss 
w tr 


x,y 

X J 

a 

r 

r 
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e 

1 

a, e, if/ 
I 

p 

T 


surface 

temperature 

temperature boundary condition 
numerical approximation of the temperature 

numerical approximation of the temperature defined by equations (3.1 la) or (4.10a) 

initial temperature 

perturbation temperature 

vector matrix containing nodal temperature 

reference temperature 

time 

final time 
initial time 

dummy integration variable 
velocity in index notation 
volume 

domain weight function 

boundary weight functions 

fundamental solution for steady-state conduction 

fundamental solution for transient conduction 

spatial coordinates 

spatial coordinates in index notation 

thermal diffusivity 

domain boundary 

integration limit 

Dirac delta function 

error residual 

transformed temporal coordinate 
polar transform variables 
transformed spatial coordinate 
solid density 
time 
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*jk 

viscous dissipation in index notation 

<Pe 

solution variable interpolation function 

Xi 

local transformed temporal coordinate 

Q 

domain 

CO 

parameter, defined by equation (3.33) 

Subscripts 


e 

denotes element 

i, j, k, m 

indices 

P 

denotes source point 

q 

denotes node point 

r 

denotes boundary 

Q 

denotes domain 

Superscripts 


ss 

denotes steady-state 

tr 

denotes transient 




This paper presents BEM algorithms for two-dimensional, steady-state and transient, heat 
conduction. The algorithms are specifically designed for efficient coupling with CFD. The presen- 
tation includes a brief tutorial on the BEM for those unfamiliar with the technique. Then the details 
of the current algorithms are presented. The mathematics are overly detailed because the paper is 
intended as a technical reference manual to document the algorithms. The document is organized as 
follows: The governing equation of heat conduction is given in Chapter 2, the steady-state algorithm 
with numerical examples in Chapter 3, and the transient algorithm with numerical examples in 
Chapter 4. 


2 Governing Equations 


Consider a material volume in space with no internal heat generation. Applying the conservation 
of energy, the time rate of change of the volumetric internal energy is equal to the net flux of energy 
through the bounding surface; in index notation, this relation is 

J e t dV + j" [e t uj + +qj jnj dS = 0 (2.1) 

Volume Surface 

where e t is the total energy per unit volume, uj is the velocity of the particles, is the viscous 
dissipation, qj is the heat flux, n j is the outward unit normal, and t is time. In a solid, the material 
velocity is zero; therefore, the convection and viscous dissipation terms are zero. Furthermore, 
since the kinetic energy is zero, the total energy comprises only the internal mode. Applying the 
Divergence Theorem, equation (2.1) becomes 

2 f edV+ f ^-dV = 0 (2.2) 

dt J J dxj 

Volume Volume 

where is the internal energy and xj are the independent spatial variables. Expressions for the 
internal energy and the heat flux are now established. The internal energy of a solid is empirically 
determined through the heat capacity, denoted as c, defined as the change in heat content of the solid 
per change in degree of temperature. In general, c is determined experimentally over a range of 
temperatures. The relation is 


_ 1 de 
C = ~p~dT 


(2.3) 


The internal energy is found by integrating cdT from a given reference temperature, T 0 , to the 
temperature of interest. 


T 

e = p Jc dT (2.4) 

To 


3 



If c is independent of temperature, then the internal energy simplifies to 


(2.5) 


e - pcT 


In general, the heat flux comprises radiative and conductive terms. Neglecting radiation, the heat 
flux is given by Fourier’s law, which states that the conductive flux is proportional to the tempera- 
ture gradient. In index notation, Fourier’s law is 


% 


= -K 


jk 


ST 

dxj 


( 2 . 6 ) 


where K k is the thermal conductivity. In general, the thermal conductivity is a tensor quantity that 
is a function of position, direction, and temperature. For this analysis, however, the thermal conduc- 
tivity is considered to be a scalar quantity denoted as k s . Thus, modifying equation (2.6) for scalar 
conductivity and substituting equations (2.5) and (2.6), the energy conservation law (eq. (2.2)), 
becomes 


I 


Volume 


dT d 

pc T + ^ 


( 


-k 




dr_ 

S dXjj 


dV = 0 


(2.7) 


In equation (2.7), the time derivative is pulled inside the volume integral; this operation is valid as 
long as the limits of integration are time independent, i.e., the domain is fixed. Letting the volume 
shrink to zero, the differential form of the heat conduction equation is derived as 


pc i :=v * ( * iV7 ' ) 


( 2 . 8 ) 


For a constant thermal conductivity equation (2.8) reduces to 


— = aV 2 T 
dt 


(2.9) 


where a = k s /(pc) is the thermal diffusivity. At the steady state, the time derivative of temperature 
is zero and the heat conduction equation reduces to Laplace s equation, given as 

v 2 r = o (2.10) 


In summary, the transient heat conduction equation given by equation (2.8) is valid for a solid 
material with a constant c. Further simplification of a constant scalar conductivity yields equa- 
tion (2.9), which is the transient conduction equation expressed in terms of the thermal diffusivity. 
Finally, in the steady state, equation (2.9) reduces to Laplace’s equation, given by equation (2.10). 
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3 Steady-State Boundary Element Algorithm 


This section develops the boundary element procedure for solving steady-state heat conduction. 
For completeness, a derivation of the boundary integral equation is presented. The presentation is 
self-contained, yet is only cursory in detail. Definitive presentations on the BEM are found in 
references 4 through 6. Regardless of the potential insufficiencies, the derivation of the boundary 
integral equation is presented through weighted residual analysis. Next, the weight function, which 
appears in the boundary integral equation as the kernel of an integral transform, is chosen to be 
Green’s free space solution to the governing equation. The free space solution and its directional 
derivative are specified in this section. Then, shape functions are presented for a linear distribution 
of the dependent variables over linear boundary elements with offset nodes. Analytic expressions for 
the integral transforms are given for both singular and nonsingular integrands. Finally, test cases are 
presented. 


3.1 The Boundary Integral Equation 

The core of the boundary element method is the boundary integral equation. Equation (2. 10) is 
the starting point for the derivation of the boundary integral equation for steady-state heat conduc- 
tion. The derivation is presented through the perspective of a weighted residual analysis based upon 
presentations in references 4 and 6. Let Q be the solid domain with boundary T upon which 
equation (2.10) is valid. Furthermore, let the boundary be divided into two parts, Tj and T 2 , for 
which the following boundary conditions apply: 

r=fonr! (3.1a) 

— = qonT 2 (3.1b) 

dn 

Further subdivision of the boundary does not enhance the validity of the derivation; it only adds to 
the complexity of the algebra. 

The goal of any numerical approximation is to minimize the error in the satisfaction of the 
governing equation and the boundary conditions. Frequently, the numerical scheme is designed to 
satisfy the boundary conditions exactly, and the error is minimized on the interior. For this deriva- 
tion, however, the strict enforcement of the boundary condition is relaxed; the numerical scheme is 
constructed to minimize the error over the domain and the boundary. Let T represent the numerical 
approximation to T, and let e be the error residual. Over the domain and boundary, the error is 
given by 


e a = V 2 i 

(3.2a) 

*r, ,=r-f 

(3.2b) 
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(3.2c) 



In weighted residual analysis, the error residuals are multiplied by weight functions and integrated 
over the domain and boundary to measure the global error. The weight functions can be viewed as 
error distribution functions whose choice determines the type of numerical approximation. For 
example, the method of Galerkin is obtained by choosing the weight functions from the same class 
of functions used to describe T. An instructive presentation of the weighted residual approach and 
its connection to finite difference, finite element, and least-squares techniques is given in refer- 
ence 6. The weighted residual expression based on the error residuals of equations (3.2) is given as 

j , £ £2 WJQ+ je ri WdT+ je^WdT^O (3.3) 

q n r 2 


where W is the weight function over the domain and W and W are the weight functions over the 
boundary. Substituting in the residual expressions of equations (3.2) yields 


JVv 2 fjO+ ^[T-T)wdT+ J[ 

o. r, r 2 ^ 


Y 


w dr = o 


j 


(3.4) 


The boundary integral equation is derived by manipulating the domain integral and judiciously 
choosing the weight functions. From the weighted residual perspective, the steps of the derivation 
appear prescient in their introduction; indeed, the source of the foreknowledge is the original 
formulation from reciprocity considerations. Rizzo (ref. 7) presents an interesting and informative 
historical view of the boundary integral technique wherein the derivation is far more deductive 
than that presented herein. Nevertheless, the focus is on transforming the domain integral of equa- 
tion (3.4) into a more convenient form. This transformation is accomplished by applying Green’s 
identity to the domain integral to reduce the order of the operator on the temperature field. The 
relation is given by 

jwV 2 fdQ = -jvw*VfdQ + jwVf»ndT (3.5) 

Q Q r 


Substitution of equation (3.5) for the domain integral of equation (3.4) produces the weak form 
of the residual statement, which is the basis of the finite element method. In the weak form, the 
weight function is symmetric to the numerical solution and, depending on the choice of the weight 
function, frequently gives rise to symmetric matrices. To produce the boundary integral equation, 
Green’s identity is applied a second time to transfer the Laplacian operator to the weight function. 
Transposing the role of f and W, equation (3.5) is rearranged to yield 

Jvw»Vf^n = -jf V 2 WdQ + jf VW*n dT (3.6) 

q Q. r 
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Successive substitution of equations (3.6) and (3.5) into equation (3.4) yields 


jfV 2 WdQ- 

Q. 


jfVW»ndT + ^WVf*ndr+j(i-T ) )Wdr + 

r r G 


[\df _ 

J A"* 


A 


\w dr = o 


J 


(3.7) 


By consecutive application of Green’s identity, the Laplacian operator in the domain integral has 
been transferred from the numerical solution to the weight function. This formulation is termed the 
inverse problem , and upon its derivation, attention is turned to the weight functions. 

Up to this point, the only limiting assumption is that the weight function, W, must be twice 
differentiable in order to apply Green’s theorem consecutively. There is complete freedom in 
selecting the boundary weights; they are chosen such that 

W =^on Tj (3.8a) 

an 

W=-WonT 2 (3.8b) 

Substituting equations (3.8) into equation (3.7) and noting the cancellations in the boundary integrals 
results in 

\TVW*hdV+\wS/T*hdr-\T^dr+\qWdr = 0 (3.9) 

q r 2 G G r 2 

By consolidating the notation, the boundary integrals can be simplified to yield 

dr+\w^dr = o (3.io) 

a r r 


where 


J T over Q and on T 2 
[T on Tj 


df 

dn 


dT 

\ dn 


on Tj 


on r 2 


(3.11a) 


(3.11b) 


It is instructive to review the steps leading to equation (3.10). First, a weighted residual statement 
is written with separate weight functions for the boundaries and the domain. Green s theorem is 
applied to yield the inverse problem shown in equation (3.7). Then the boundary weight functions 
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are selected in terms of the interior weight function to produce cancellation in the boundary 
integrals. Finally, a judicious variable change simplifies the integral equation. 

The remaining task is the selection of the weight function; a profitable choice is Green’s free 
space solution to the governing equation. Green’s function is a fundamental solution to the govern- 
ing equation subject to a unit impulse forcing function. The fundamental solution for steady-state 
conduction, denoted by W ss , satisfies 


V 2 W ss =d{q-p) (3.12) 

where S is the Dirac delta function, q is a position vector to any point in the domain, and p is the 
position vector to the point source. The precise mathematical definition of the Dirac function is 
ambiguous, but its critical property is that, for a function F(x), the integral of the product of 
F(x)S(x- xq) satisfies 

Jf(x)<5(x-xo)^Q = F(xo) (3.13) 

Q 

In some presentations, the Dirac function is defined by equation (3.13), and in some instances, 
equation (3. 13) is a property of the definition; a more extensive discussion of the Dirac function is 
given in reference 8. In any case, the operation of the Dirac function in lieu of V W ss in the domain 
integral isolates the value of the temperature at the source point. The domain integral is effectively 
eliminated. The boundary integral equation becomes 

CJ D - dT+ dT = 0 (3.14) 

r r 

where f p is the numerical approximation of the temperature at the source point p. The coefficient 
C p is a function of the included angle exposed to the interior at the source point. Details of the 
derivation of C p are given in reference 4. 

The boundary integral equation is the core of the boundary element method. By choosing the 
fundamental solution as the weight function, domain integration has been eliminated; observe from 
equation (3.14) that only boundary integrals appear. The result of this development is a numerical 
procedure where the nodal points are located only on the boundary. The boundary-only character 
stands in contrast to finite difference or finite element techniques that require a complete domain 
discretization; the benefit is a reduction in dimensionality. Furthermore, with the BEM any subset 
of the interior solution can be calculated to any desired resolution based on the computed boundary 
solution. 

The general outline of the BEM is as follows: The boundary is discretized into elements that can 
be of any shape, but typically are polynomials as in finite element procedures. The solution variables 
are assigned an interpolation function based on nodal points distributed over the element. The inter- 
polation function defines the distribution of the solution variable over the element. Originating at the 
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source point, the fundamental solution and its directional derivative are kernels of an integral 
transform of the prescribed interpolation functions. After dividing the boundary into elements, the 
boundary integral equation becomes 

v,-ijn^+xj w ”§ < 3 - 15) 

e r- e p 

I e l e 

where e denotes an individual element, so f e and dtjdn denote the distribution of the dependent 
variables over element e. Typically, the integrals are computed numerically using Gaussian 
quadrature. 

Equation (3. 15) is written for each node to form a system of linear equations. The system is 
expressed in matrix notation as 

H jt T,+G; t Qj=0 (3.16) 


where T and Q j are vectors containing the nodal temperatures and temperature gradients, 
respectively, and^ H jk and G jk are matrices containing the influence coefficients resulting from 
the integral transform. After segregating the known and unknown dependent variables based on the 
boundary conditions, the linear system is solved to yield the complete solution on the boundary. 

The algorithms presented in this paper employ Gaussian elimination with partial pivoting for direct 
inversion of the system matrix. To compute the interior solution, equation (3.15) is applied with the 
source point located at the interior point of interest. Since the solution on the boundary is completely 
known, the boundary integrals are computed directly without a matrix inversion. 

The remaining ingredients of the numerical recipe are the definition of the fundamental solution 
and its derivative, the geometrical definition of the boundary element, and the prescription of the 
approximating interpolation function to the unknown solution variables f and df/dn. The combi- 
nation of these ingredients differentiates particular boundary element algorithms. The algorithms 
presented here are specifically designed for efficient coupling with a CFD flow solver. Serendipi- 
tously, an analytic solution of the integral transforms is achieved with the chosen combination of 
ingredients. 


3.2 Fundamental Solution for Two-Dimensional, Steady-State Heat Conduction 

The fundamental solution for the two-dimensional Laplace’s equation is given by 

W ss = — In — 

2 71 r 

where r is the distance from the source point to a point in the domain; it is expressed as 

r = ||r|| = |g-p|| 


(3.17) 


(3.18) 
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The directional derivative is given as 


dW ss = 1 r*n 
dn 2 k r 2 


(3.19) 


Both the fundamental solution and its derivative are singular at r = 0, so care must be taken when 
integrating near or through the source point. 

3.3 Element Shape Function 

The next ingredient to the numerical recipe is the definition of the boundary element shape. The 
element shape is distinct from and prescribed independently of the dependent variable distribution. 
The element shape is a geometrical attribute that determines the distance function, the outward 
normal, and the integration path of the contour integrals. On the other hand, the dependent variable 
interpolation function defines the distribution of the dependent variables over the element. The two 
functions are constrained differently; the element shape is determined by the physical domain, 
whereas the interpolation function is governed by the variation of the solution over the boundary. 
The two functions combine to determine the accuracy of the algorithm. The element shape function 
is given in this section while the variable interpolation function is described in the next. 

This boundary element algorithm is specifically designed for coupling with CFD codes that 
employ finite difference or finite volume techniques. For both structured and unstructured grids, 
these CFD techniques assume linear segments between grid points. Thus, linear boundary elements 
are selected in order to ensure one-to-one correspondence of the boundary element grid to the CFD 
grid. With linear elements, conservation is easily satisfied since interpolation is not required to mate 
the two domains. Furthermore, employing linear boundary elements creates two simplifications. 
First, the distance function, r, is prescribed analytically between any arbitrary point and line seg- 
ment. Second, the outward unit normal to a linear element remains constant; consequently, r • h , 
which arises from the directional derivative of the fundamental solution, is pulled out of the 
integration. These simplifications allow analytic solution of the integral transforms. 

The component expressions for r and r •ft for a linear element are now presented. In two- 
dimensional space, the source point can be oriented with respect to a linear segment in one of three 
possible ways. Each orientation results in a different analytic expression and is addressed separately. 
In Orientation 1, the source point is not collinear with the line segment; thus, r • h * 0. For Orien- 
tation 2, the source point and element are collinear but the source point does not lie on the segment 
itself. Finally, in Orientation 3, the source point lies on the element; thus, the element contains an 
integrable singularity. In both Orientations 2 and 3 , r • n = 0. 

3.3.1 Orientation 1 Distance Function- Figure 1 displays the notation and schematic of a source 
point and linear boundary element in Orientation 1; r is the position vector from the source point to 
a point on the element; r 12 is the position vector from the first point of the line segment to the end 
point; and r pi is the position vector from the source point to the first point of the line segment. The 
element is mapped into a linear segment of unit length through the transformation 
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(* 2 > > ! 2 ) 



(*/>. ») 


\ 


Figure 1. Source point and linear boundary element segment in Orientation 1. 


x = x\ +£(*2 -*i) (3.20a) 

}’ = Ji + £(J2 ~Jl) (3.20b) 

where (xj,jj) and are the start and end coordinates, respectively, of the line segment 

defining the element, and E, is the transform variable with range 0 < E, < 1 . The square of the 
distance function varies as a quadratic function of the transform variable according to 


r~ = r • r = A + BE,+ C£ 2 

(3.21) 

where 


II 

■£*' 

• 

(3.22a) 

B = 2r pl »r n 

(3.22b) 

C = r n *r | 2 

(3.22c) 

The coefficients A, B, and C are solely functions of the geometry; they attain different values for 


each element and source point combination. 
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The magnitude of the dot product of the position vector and the element outward unit normal can 
be expressed in terms of the distance function coefficients as 


r • n = 


2s[C 


(3.23) 


where 


Q = 4AC-B 2 (3.24) 

By ordering the elements in a counterclockwise fashion by convention, the sign of the dot product is 
determined by 


sign(r • h) = sign(r pl x r 12 ) 


(3.25) 


Note that the dot product becomes zero when Q is zero. 


3.3.2 Orientation 2 Distance Function- In Orientation 2, the source point is collinear with the 
element. As will be shown later, the analytic expressions of the integration become undefined when 
the source point and line segment are collinear, i.e., when Q — 0. This degeneracy is avoided by 
redefining the distance function. Since the source point and element are collinear, the location of the 
source point can be expressed in terms of the element transform variable £. Let £ p be the position 
of the source point. Recall that the element is defined over the range 0 < £ < 1 ; thus, if 0<£ p <1, 
then the source point lies on the line segment resulting in Orientation 3. In Orientation 2, the source 
lies off the element, so t; p > 1 or £ p < 0 . In either case, is determined as 



(3.26) 


The distance function becomes 


(3.27 

| a/A + $4C for E, p < 0 

3.3.3 Orientation 3 Distance Function- In the third orientation, the source point lies on the 
element and hence the integral transform is singular at a point along the path of integration. Never- 
theless, the integral transform exists in the Cauchy principle value sense. In order to compute the 
principle value, it is convenient to set the origin of the transform variable at the source point. Such 
a transformation is dependent on the location of the source point with respect to the grid points; 
consequently, it is dependent on the solution variable shape function. Therefore, the discussion of 
the distance function for Orientation 3 is delayed until after the presentation of the solution variable 
interpolation function in the next section. 
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3.4 Solution Variable Interpolation Function 


The final ingredient of the numerical recipe is the dependent variable interpolation function. 
Barring an algorithm designed for a specific application, it is common to set the interpolation 
function to the same order as the element shape because, in general, the least accurate function will 
constrain the global accuracy. Thus, the distribution of the solution variables is selected to be linear. 
Two node points are required on each element to properly define a linear distribution. Normally in 
finite difference and finite element methods the node points are coincident with the grid points or, 
in other words, the node points are located at the ends of the linear element. For piecewise linear 
segments, however, the element normal direction will most likely change discontinuously from one 
element to the next. Consequently, when the nodal points are coincident with the grid points, the 
node normal direction, which is used to define dT/dn, becomes nonunique at the node. There are 
procedures to account for this effect, which amount to substituting an auxiliary equation when 
relating the upstream and downstream temperature gradients of a particular grid-node point (for 
example, see ref. 9). In the present work, the node points are variably offset from the ends of the line 
segment. The placement of the node away from the grid point uniquely defines the element normal 
at the node. The drawback of the offsetting strategy is an increase in the number of nodes since, in 
general, two unique nodes per element are required, whereas a node located at the grid point is 
common to the two adjoining elements. This offsetting procedure is not required in all circum- 
stances. For instance, if two adjoining elements are collinear, then the node point is placed at the 
grid point without concern. Since minimizing the number of nodes to reduce the computational 
effort is desirable, a shape function that accounts for variable positioning of the node is described. 

Figure 2 shows an element and the position of the nodal points. The nodes are offset from the 
beginning and end grid points by transformed distances a and b, respectively. If the value of a or b 
is zero, then the grid and node point are coincident. Let represent either of the solution variables 
f or df/dn over element e. Also, let 0 e l and 2 be the values of the solution variable at the 
beginning and end nodes of element e, respectively. The solution variable is given by the linear 
function 


h = — ‘ ■ ^ [fe,l( 1 -fr-{)+fe.2«-«>] (3-28) 

Although not indicated by a subscript e, the values of a and b vary from element to element, 
depending on the local geometry and boundary conditions. 

• Grid Point 
o Node Point 

(xi.yi) *'•' (*2,ys) 

> 0 ■ O • 

II II 

£ = 0 Z = a $ = \-b £ = 1 

Figure 2. Linear boundary element definition. 
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3.5 Integral Transforms 


Definitions of the weight function, the element shape, and the variable distribution function are 
now complete. These elements are combined to derive analytic expressions for the integral trans- 
forms. Two integrals are computed for each element and source point combination. The integrals 
are mapped into the transform variable £ and are given as 



dW 


SS 


dn 


-dr = 


f \ySsdTe_ 

J dn 


dr = 



(3.29a) 


(3.29b) 


where l e is the length of the element. Recall that the distance function and dot product contained in 
the fundamental solution and the solution variables have all been previously expressed in terms of 
£ . The integral expressions are analytically derived for each of the three possible orientations. 

3.5.1 Orientation 1 Integrals- In this orientation, the fundamental solution, with the expression 
for the distance function given by equation (3.21), is combined with the shape function of equa- 
tion (3.28) to produce the integral transforms. The integrals of equations (3.29) are of similar form 

i = w - «> { fgJ l (1 ~ bvr - /? K 2 1'" ~ a/ '"l } (3 ' 30a) 


X 

Jo 


i A w 


SS 


dn 


d£ = 


Vc 

An(\ -b-a) 


( *Vr 




e,2 


(3.30b) 


The terms if denote component integrals resulting from algebraic manipulation of the integrand. 
The superscript ss indicates steady-state integral components and is used to distinguish those 
integrals from similar transient integrals presented in the next chapter. The component integrals 
are derived analytically as 


j SS 

l \ 




f'dS_ J_ 
hr 2 'JO 

+ B + C| - ln|4- S/f ) 


tan 


-l 


f 


2C+JB 

4q 


tan 



(3.31a) 

(3.31b) 


r ss 

It 




In \r z 




inW+^/f-2 


(3.31c) 


14 



(3.3 Id) 


If = £ 4 ln|r 2 | dS = i ln|A + B+ C] - 1 / 5 “ - Clf 


jSS 

l 6 


= r i £^ = in 

"Jo r 2 cb 


- - 5/f - A/f 


(3.31e) 
(3.3 1 f) 


3.5.2 Orientation 2 Integrals- In the second orientation, r • n = 0 and therefore the integral of 
equation (3.29a) is zero. The expression for equation (3.29b) is given by equation (3.30b) with the 
integrals if and /" derived from the distance function given by equation (3.27). The integrals are 
given as 


If = 2 






CO 


ln| |VA + cy| - 1 - ^-ln|VA| 


! 4 = f 1 - £) + H + ^ ~ \ 


where 


J VC for £ p <0 

|-VC for ^>1 


(3.32a) 

(3.32b) 

(3.33) 


3.5.3 Orientation 3 Integrals- As in Orientation 2, r *n = 0; therefore, the integral of equation 
(3.29a) is zero. In the third orientation, the integrand of equation (3.29b) is singular at the source 
point, but the integral exists in the Cauchy principle value sense. The derivation is accomplished by 
recasting the distance and interpolation functions in terms of a coordinate system originating at the 
source point. It is possible for the source point to be located at either of the two nodes on the 
element, and the singular integration is similar for both possibilities. A slight change in notation 
condenses the equations of the two possibilities into a single expression. Let subscript p denote the 
value of the gradient at the source node and let subscript q denote the other node. In this notation, 
a and b are the distances separating points p and q from the grid point, respectively. This notation 
departs from the previous definition, which associated a and b with the ordering of the element 
nodes. In this new definition, the ordering is arbitrary and the equations are applicable for either 
node point as source. The integral is given as 


J 


w ss^e. dr = je_ 
dn 2 n 


3f„ dt a 

-f -(Aj— A 2 ) + -f-A 2 
dn dn 


(3.34) 
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where 


A\ =ln|/ e |-l + (l-a)ln|l-a| + aln|a| 


(3.35a) 


Aj — ■ 


1 


2(1 — b — a) 


(1 - 2a)^ln|/ e | - + (! - a) 2 ln|l - a\ + a 2 ln|a| (3.35b) 


When a = 0, the previous expressions for the coefficients A\ and Ai are undefined because of the 
natural log function; however, the limiting form is found to be 


A 1=141-! 


(3.36a) 


A~> — ■ 


1 


2(1 ~b) 




(3.36b) 


The numerical recipe is complete. Numerical solutions are generated by looping over all node 
points and elements to determine the orientation. Then the appropriate analytic solutions are used to 
compute the elements of the global system matrix. Subsequent inversion of the system matrix 
produces the boundary solution. 


3.6 Numerical Results 

The previously developed steady-state algorithm is tested by comparing the numerical solution 
to three separate analytic solutions. The first test case is a cylindrical geometry with imposed surface 
temperatures. The second test case is a pin-cushion-like geometry formed by the intersection of four 
circular arcs. The third test case is a seven-sided nonconvex polygon. 

3.6.1 Cylinder- The first test problem is the computation of the temperature distribution in a 
cylinder with an inner radius of 1 unit and an outer radius of 2 units, as given by Becker (ref. 5). The 
temperatures on the inner and outer radii are 10 and 6, respectively. The exact solution is given by 

7 = -5.771 lnr + 10 (3.37) 

The boundary is discretized using 24 elements on the quarter plane: 8 on each of the inner and outer 
radii and 4 on each of the remaining two sides. The temperature is prescribed on the inner and outer 
radii and symmetry boundary conditions ( df/dn = 0) are imposed on the two sides. The steady- 
state temperature contours of the boundary element solution are compared to the analytic solution 
in figure 3. The contours demonstrate the radial symmetry of the numerical solution both on the 
interior and along the symmetry boundary. Furthermore, it is seen that the BEM solution compares 
very well to the exact solution. The jaggedness of the inner and outer radii curves on the boundary 
element side of the plot is a result of the piecewise linear segments used to describe the circle. 
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Figure 3. Comparison of analytic and boundary element steady-state temperature contours in a 
cylinder. 


3.6.2 Pin Cushion- The second test case is found in reference 9; it is a pin-cushion geometry 
defined by the intersection of four circular arcs centered at (±1, ±3) with a radius of 3.64, as shown 
in figure 4. The geometry contains four comers located at the intersection of the circular arcs. 
Potential deficiencies of the offset node strategy are most likely to appear at such surface contour 
discontinuities. The exact solution is taken to be 


T(x,y) 


x 2 - y 2 + In 


[(x-2) 2 +(y-2) 2 ] 2 


(3.38) 


Contour lines of the exact solution are shown in figure 4. It should be noted that the temperature 
field is not symmetric about any axis, consequently, the entire domain must be computed. 

The boundary is discretized using 8 elements of equal length on each side for a total of 
32 elements. The analytic temperature is imposed on the boundary. The computed surface tempera- 
ture gradient is compared to the analytic solution in table 1 . Comparisons are made upstream and 
downstream (moving counterclockwise) at the four comers of the pin cushion; typically, the errors 
in the computed gradients are greatest at the comers. The BEM gradient is within 5 percent of the 
analytic flux at the selected points despite the potential inaccuracies generated by offset nodes. A 
similar comparison by Kassab and Nordlund (ref. 9) shows their method to produce errors less than 
0. 1 percent. The improved accuracy of their method results from the use of quadratic elements, 
which not only model the surface properties to higher order accuracy but are also able to reproduce 
the exact surface geometry ; the linear elements used in this study generate error from approximating 
the circular arc shape as a series of line segments. 

For the solution reported in table 1, the nodes are offset from the end points by 10 percent of the 
element length (a = b = 0.1). The solution is fairly independent of the offset distance in the range 
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Point 

Analytic dT/dn 

BEM dT/dn 

Comer 1 upstream 

1.1959 

1.1758 

Comer 1 downstream 

0.9494 

0.9094 

Comer 2 upstream 

-1.2495 

-1.2489 

Comer 2 downstream 

-1.1423 

-1.1280 

Comer 3 upstream 

0.9735 

0.9501 

Comer 3 downstream 

1.6489 

1.5978 

Comer 4 upstream 

-0.6427 

-0.6370 

Comer 4 downstream 

-0.8186 

-0.8151 


lErrorl 

0.0201 

0.0400 

0.0006 

0.0143 

0.0234 

0.0511 

0.0057 

0.0035 


Percent error 


-1 

-4 

-0 

-1 

-2 

-3 

-0 

-0 


0.01 < a,b < 0.25. In general, experience has shown that a = b = 0.1 is the most reliable offset value; 
all the results reported herein are computed using a 10-percent nodal offset. Also, the integral trans- 
forms are computed using 16-point Gaussian quadrature to verify the exact integration equations 
outlined previously. The gradient values generated using Gaussian quadrature prove to be within 
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4 significant digits of the exact integration results. Additionally, the interior temperature field is 
computed based on the surface solution. Temperature contours compare within the plotting accuracy 
of figure 4. 

The spatial accuracy of the algorithm is measured by computing an error norm for successive 
increases in the number of elements. The L 2 error norm is shown in figure 5 as a function of the 
number of boundary elements for both the boundary gradient and the interior temperature. The error 
norms of figure 5 are normalized by the respective values of the error using 16 elements. The results 
show the calculation of the boundary gradient to be first-order accurate while the computation of the 
interior temperature is second-order accurate. This result is consistent with the accuracy of linear 
elements employed in the finite element method. 



Figure 5. Computed error of pin-cushion solution. 


3.6.3 Nonconvex Polygon- The final case is a nonconvex polygon, also extracted from reference 9. 
The geometry, which is shown in figure 6, is designed to test a range of comer angles. Also shown 
in figure 6 are contours of the exact solution, given by 

T(x,y) = sin(x)cosh(y) (3.39) 

The boundary is discretized using 26 elements matching that of reference 9. The analytic tempera- 
ture is imposed on the boundary. A comparison of the computed boundary temperature gradient is 
given in figure 7. The abscissa of the plot is surface length starting at ( x,y ) = (1,1) and proceeding 
counterclockwise around the domain. Also displayed are the grid point locations. The discontinuities 
in the gradient values result from the discontinuous change in the surface normal between sides of 
the polygon. As seen in figure 7, the computed temperature gradient is within 5 percent of the 
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analytic value on the boundary except at the comer points. The worst comparison is at the 90° comer 
at (x,y) = (2,1). It is unclear why the maximum error occurs at this location. The solution is not 
particularly ill-behaved at the comer, suggesting that the effect is geometry governed. It is believed 
that this result is particular to the offset node strategy. In reference 9, the error at the comers corre- 
lated to the turn angle with the 1 17° comer at (x,y) = (2,1.5) showing the greatest error; no such 
correlation is found in the solution shown in figure 7. In any case, the efficacy of the solution 
procedure is demonstrated. Finally, although not shown but worthy of mention, contours of the 
interior temperature field compare to the analytic solution within the plotting accuracy shown in 
figure 6. 


4 Transient Boundary Element Algorithm 


The development of the boundary element method for transient heat conduction follows the same 
path as the steady-state presentation given in Chapter 3. In principle, the procedure is the same as the 
steady-state, however, the mathematics are ubiquitous, and the general procedure is easily obscured 
by the mathematical details. Referral to Chapter 3 may illuminate Chapter 4. 

4.1 The Boundary Integral Equation 

The derivation of the boundary integral equation commences with the transient heat conduction 
equation expressed in terms of the thermal diffusivity given by equation (2.9) and repeated as 
follows: 


ST 

dt 


with time-varying boundary conditions given by 


dT 

dn 


= 0 in Q 

(4.1) 

on Tj 

(4.2a) 

onr 2 

(4.2b) 


The transient temperature field is split into an initial temperature field, 7), and a perturbation 
temperature, T ' . The relationship expressed in terms of functional dependence is 

r(x ; ,o=n*j.O +*>(*/) (4.3) 

Furthermore, it is assumed that the transient calculation is computed from a steady-state initial 
condition, i.e.; V 2 7) =0. Accounting for these relationships, equation (4.1) becomes 

— -«V 2 r = 0 in £2 (4.4) 

dt 
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Equation (4.4) is the basis for the development of the boundary integral equation. Since the 
governing equation is written in terms of the perturbation temperature, boundary and initial condi- 
tions must be expressed in like terms. The initial conditions are T ( xj,t = 0) = 0 by definition. The 
boundary conditions are recast in terms of the perturbation temperature as 


T' = T(t)-T, onT x 

dT' . dr, „ 

a— = q(t)-a—±- onT 2 
dn on 

Next, the domain and boundary error residuals are constructed: 

dT „2-f 

e Q =-^~aV T 


fiT, =T-T + T, 


dr _ dr, 

er2= «__, +a _ 


(4.5a) 

(4.5b) 


(4.6a) 

(4.6b) 

(4.6c) 


where T is the numerical approximation to V . Integrating over the domain, boundaries, and time, 
the weighted residual statement becomes 

Jjf^_aV 2 fy^a^T + J J(f-f + 7))iyjrrfT + J jfa^-9+a^- W dT d? = 0 (4.7) 

rD 2 tTi tT 2 V 


where W is the interior weight function, W and W are the weight functions on the boundaries, and 
r represents the time domain. Up to this point, the procedure is analogous to the steady-state deriva- 
tion with the additional dimension of time. The inverse problem is derived by using Green’ s theorem 
to transfer the Laplacian operator from the temperature field to the weight function. In addition, the 
integral of the time derivative of temperature is transformed through the relation 



— WdQdt = 

dt 


tQ 




f 

J 'o 


' f Ar 

dt 


dn 


(4.8) 


where r 0 and tf are the initial and final time, respectively. The boundary weight functions are 
selected in terms of the interior weight function as 


u7 dw „ 

W = -a— —on f, 

dn 

(4.9a) 

W = W on r 2 

(4.9b) 
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Furthermore, the notation is consolidated through the variable change given by 


Y _\r over ^ and on T 2 
J ~\ T-Tj on Tj 


dt 


31 

a— on l] 

an 

- dT, _ 
q - a— L on T 2 
on 


(4.10a) 


(4.10b) 


After the application of Green’s theorem, the substitutions of equations (4.8 through 4.10), and some 
algebra, the weighted residual statement becomes 


f TW 

— TW 

\ 

dQ- f [ f(^ + aV 2 W 

r\ v 


'o, 

J J V dt J 

ST 


dQd t 


a j{ w ^ drdr + a jj^^ drdT = 

r r x r 


0 


(4.11) 


In order to remove the second domain integral, the weight function for the transient equation, 
denoted as W tr , is selected to satisfy the following: 


dW tr 

dt 


+ aV 2 W tr = -8(r)S(t) 


(4.12) 


The right-hand side parameters of equation (4.12) are Dirac functions,' which represent a unit 
impulse forcing function in space and time. In addition, Wrobel and Brebbia (ref. 10) have shown 
that 


J 


TW 


tr 


dQ = 0 


(4.13) 


Q '/ 

Applying equations (4.12 and 4.13) to equation (4.1 1) yields the transient boundary integral 
equation, given by 


~[fW tr dQ. + C p t p -a^W tr ^-drdx + a^f^^-drdx = Q) (4.14) 

a r ° rr ir 


Other than the domain integral at t = t 0 , the transient boundary integral equation is identical in form 
to the steady-state boundary integral equation (3.14). The differences are the weight function, the 
additional dimension of time, and the presence of the domain integral, which appears to prevent 
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boundary-only discretization. Nevertheless, the numerical methodology is identical to the steady- 
state algorithm. Moreover, it is shown in Section 4.3 that the calculation of the domain integral can 
be circumvented to recover the boundary-only character of the algorithm. Before that analysis, the 
transient fundamental solution is introduced. 

4.2 Fundamental Solution for Two-Dimensional, Transient Heat Conduction 


The fundamental solution satisfying equation (4.12) is given by 


W tr = 


Ana{tf -t) 


exp 


— r 


.2 


4 a(tf-t) 


H(t f -t ) 


(4.15) 


where H(tf - 1 ) is the Heaviside function. The fundamental solution is singular when r = 0 and 
t = tf . The directional derivative is given by 


3W tr 

dn 


-(r • n) 

2 2 eX P| 

8 KGC Z {tf-t) 1 y 


-r 


4 




(4.16) 


Both the fundamental solution and its directional derivative are functions of r and ( tf-t ); i.e., 
W tr = W tr {r,tf - 1 ) . This functional dependence is important in the development of the time-step 
procedure presented in the next section. 


When computing the transient integrals, it is convenient to transform the fundamental solution to 
normalized coordinates. The spatial coordinate is mapped to the parameter % previously introduced 
in Chapter 3. Recall that the distance function, r, is expressed in terms of £ . Time is transformed 
into the coordinate 77 , which originates at the source time level, tf, and is scaled by the time 
interval ( tf - 1^). The transformation is backward in time and is given by 


77 = 


IfzL 

tf-t 0 


(4.17) 


The transform variable 77 ranges from 0 at t = t f to 1 at t = 7 0 . The fundamental solution and its 
derivative are written in terms of the transformed variables as 


W tr (Z,T],tf-t 0 ) = 


1 -rhb 

exp 

4na(t f - 7 0 )77 \^4a(tf-t 0 )ri 


H(v(t f -t 0 )) 


(4.18) 


dW°\^ltf -7 0 ) _ -(r»«) 


O — T®Xp 

dn Stta (tf -torn { 


-r 2 (£) 

4a(t f -t 0 )ri 


HWtf-to)) 


(4.19) 
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The fundamental solution and its derivative are functionally dependent on £ , r \ , and the time 
interval (tf - to) from the temporal coordinate transformation. 

4.3 Time-Step Procedure 

Consider the boundary integral equation expressed over the temporal domain from /q to tf, 
given as 


- \fw tr (rj f -t ) 

Q t 0 


dQ+C p T p 




~dW tr (r,tf -t) 
dn 


dT d T = 0 


(4.20) 


The functional dependency of the fundamental solution is included in the previous expression to 
accentuate the source point, which is always located at tf. After the boundary is discretized, the 
integral equation can be written in terms of coefficient matrices of the temperature and gradient 
vectors. The form of the equation is dependent on the order of the temporal accuracy. For discussion 
purposes, the temporal accuracy is assumed to be 0(1) in order to simplify the equations. In practice, 
the algorithms presented in this paper employ a linear interpolation function, the details of which are 
presented in Sections 4.4 and 4.5. As before in the steady-state analysis, the notation is switched to 
matrix form for compactness. Equation (4.20) becomes 


Dj (Tj )* r / - { 0 ) + Hj* (tf - 1 0 )T 'j (tf) + Gj k (tf -to)Qj(tf)-0 (4.21) 


where Hy and Gj k are the coefficient matrices, identical in function to those of equation (3.16), 
and Dj is a vector containing the domain integral, which is a function of the temperature field at 
tQ and the time interval (tf - to)- The parenthetical expressions after Tj and Qj indicate time level, 
while those after Dy, H j k , and Gy* denote functional dependency. As will be shown later, the 
temporal portion of the integral transforms admit analytic solutions that are functions of the time 
interval (tf — to)- Thus, the coefficient matrices Dy, Hy* , and G y* are denoted as functionally 
dependent on the time interval. Although not explicitly indicated, the matrices are also functions 
of the geometry. 

Wrobel and Brebbia (ref. 10) present two methods for advancing the boundary integral equation 
in time. In the first approach, the time integration is initiated from the previous time step; this pro- 
cess requires domain integration and, consequently, a domain grid. In the second approach, domain 
integration is avoided by initiating the integration for every time step from r 0 - The two approaches 
compel different storage and computation strategies. Both methods are now discussed. 

Let the time domain be divided into intervals with the time at the end of each interval denoted by 
t m . In the first approach, the solution from t m _ ] to is given by 
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Dy(Ty(f m _i),? m _ ^m-l) + Hyjt(^ WJ ~ ^m-l)Ty( f m) + G^(? m ^m-l)Q;'(*m) 0 (4-22) 

Equation (4.22) is essentially a rewrite of equation (4.21) with a change in the time interval from 
(tf -to) to (t m — t m _i). To advance the solution in time, the coefficients in H j k and G jk are 
computed based on the local time step, (t m - t m -\)- The domain integral is computed based on the 
temperature field at the previous time step. For varying values of the time step, H j k and G jk need 
not be retained. If the time step is constant, however, the coefficient matrices are a function of a 
fixed geometry and a fixed time step; therefore, they can be computed once and stored. This 
procedure is very efficient except for the requisite domain integration. 

In the second method, domain integration is avoided by writing equation (4.21) from t 0 to t m . 

At the initial time, T,(/ 0 ) = 0 by definition (recall that T ) contains the perturbation temperature and 
not the absolute temperature, which can be nonzero), consequently, the domain integral is zero. 
Equation (4.20) becomes 

m m 

X H jk {t m - tf_i)T ; (r,) + XG j k {t m - ti-OQjiU) = o (4.23) 

(=i (=i 

The summation appearing in equation (4.23) results from dividing the time integral into intervals 
corresponding to the time-step distribution. The integrals are computed using the transient history 
of the dependent variables. The coefficient matrices must be calculated for every combination of 
( t m - tj) for every time step. For example, assume the solution has been computed up to time t \ . 

To advance the solution to time t 2 , coefficient matrices must computed for time intervals of (t 2 - to ) 
and {t 2 - r t ). The terms of the first summation of equation (4.23) are 

X H jk (h-ti-i )T j (tj ) =Kj k (t 2 -to )Tj ( t \) + n jk (t 2 - h )Tj (t 2 ) (4-24) 

i = 1 

At the next time step, the matrices must be computed for time intervals of (r 3 -t Q ), {t 2 -t x ), and 
( r 3 - 1 2 ). The terms of the first summation are 

X (^3 - ti-i )Ty (?/ ) = H j k (t 3 - 1 0 )Tj (q ) + H ;Jt (t 3 - q )T / - (t 2 ) + Hj k (t 3 - 1 2 )T ; - (t 3 ) (4.25) 

(=i 

Thus, for a calculation of m max nonuniform time steps, ^ + 1)/2 coefficient matrices 

must be computed to cover the possible combinations of ( t m — tj) for 0 < i < m and 1 < m < w rnax • 

In general, the value of ( t m - tj) need not repeat itself, therefore, H jk {t m - tj) and G jk (t m - tj) are 
never repeated. Storage is minimal because only one matrix at a time is needed, but the amount of 
computation is high. 

Matters are simplified somewhat for a constant time step because the number of matrix calcula- 
tions is reduced to . Furthermore, the matrices, which are a function of integer multiples of the 
time step, are computed once and stored. For a constant time step, equation (4.23) becomes 
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(4.26) 


m m 

X H jk((m - i + l)AOT;(/ f ) + XG;*((m - 1 + l)A/)Q ; (/ ; ) = 0 
i=l i=l 

Unfortunately, because the origin of the source point is always at t m , the matrix vector multiplica- 
tions of Hj k ((m - i + V)At)Tj(t l ) and Gj k ((m - i + l)At)Qj(tj) are never repeated. For example, at 
the second time step, H j k (2At) is multiplied by T^q), while at the next time step it is multiplied 
by Tj(t 2 ). Consequently, there is no possibility for computational and memory savings by storing 
intermediate vector-matrix multiplication results. 

Wrobel and Brebbia (ref. 1 1) recommend the first method, which requires domain discretization, 
for most applications. Although the first method has the overhead of the domain integral, the effi- 
ciency is independent of the number of time steps. In the second method, on the other hand, the 
efficiency decreases as the number of time steps increases; the efficiency of the second method 
quickly falls below that of the first as the number of time steps is increased. Yet, it is the author’s 
opinion that the primary benefit of using the BEM is the boundary-only discretization. If the domain 
must be discretized, perhaps the finite element method is better suited since it is conceptually less 
complicated and is amendable to increases in fidelity such as modeling nonlinearity. In light of this 
philosophical perspective, the transient algorithm presented in this paper employs the second method 
to exploit the advantage of boundary-only discretization. Of course, these brush strokes are broad, 
and certainly there may be compelling reasons to use a BEM that requires domain discretization. It 
has to be the engineer’s discretion when to use a wrench and when to use a ratchet. 

4.4 Element Shape Functions and Integral Transforms 

This section presents the numerical details of the transient boundary element algorithm. The 
transient algorithm retains the features developed for the steady state: a linear boundary element 
shape for efficient coupling with CFD, variably offset nodes for unique definition of the node 
normal, and a linear spatial distribution of the dependent variables. These components are 
incorporated into the boundary-only time-step procedure introduced in the previous section. The 
interpolation function is modified to linearly approximate the dependent variables over space and 
time. The interpolation function is inserted into the boundary integral equation to derive the integral 
expressions. The integral transforms over time and space do not admit a complete analytic solution 
except in the instance of a singularity. The integral transforms of the singular and nonsingular types 
are addressed separately. 

The boundary integral equation for the boundary-only time-step procedure is given as 
l rr r/f ir~ dW tr 

C p fp — oc jj W tr (r,t m — t)—^dT dx + a jjT-^—(r,t m -t)drdt = 0 (4.27) 

f 0 r 'of 

where t m is the time at which the solution is to be computed. It is assumed that the solution has 
been previously computed up to time t m _\ . As outlined in the previous section, the domain integral 
vanishes because the integration always commences from tQ. The time intervals are not necessarily 
constant; that assumption will be made later in the derivation. The boundary is divided into elements 
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and the time integration is divided according to the time-step distribution to yield the double 
summation in the boundary integral equation, 

-yiri ftl - Drr/fr' 

w tr (r,t m -t)^d TJr + aX J J j \n^ r (r,t m -t)drdt = 0 (4.28) 

where T l e is the distribution of f over element e from time f,_i to t h and dT l e jdn is defined 
likewise. Isolating an individual term from each of the double summations reveals the two integral 
transforms that comprise the crux of the numerical solution, 

f 1 ' f T l ZlldTdx (4.29a) 

J,_ J * dn 

* e 

f*' f W ,r (.r,t m -t)^-dTdx (4.29b) 

Jt: \ J on 

t -l— 

1 e 

The solution of the integrals is different for singular and nonsingular integrands. In the nonsingular 
case, a mixed analytical-numerical solution is prescribed. In the singular case, an analytic series 
solution is derived. 

4.4.1 Nonsingular Integrals- In this section, the equations for the integral transforms of equa- 
tions (4.29) are presented for the case of a nonsingular integrand. As will be shown, the temporal 
integration admits an analytic solution, whereas the spatial integration is performed numerically. 
First, the interpolation function of the dependent variables is developed. Let <p e represent the distri- 
bution of either f or df/dn over element e from time t { _ x to The four-point basis function is 

given as 


C p tp-af, | Ij 

i=l f . . e p 
h - 1 l e 


(1 -b-a)l 1 ( 4 . 30 ) 

+ (1 - *,)[(! - b - + (£ - a )^ 2 ('«)]} 

where Xi is a local time variable that ran S es from 1 when 1 ~ h -\ to 0 when t = t i - The parenthetical 
expressions after fe ,• indicate time level. The subscripts 1 and 2, first introduced in the steady-state 
interpolation function, indicate the two nodes of each element. Equation (4.30) is a bilinear function 
over space and time expressed in terms of the local time variable, Xi > anb tbe s P aba ^ transform 
variable £ . The interpolation function is more suitably expressed in terms of the global time 
variable T) given by equation (4. 17), which in this usage ranges from 1 at t — ?q to 0 a t t — t m . 
Equation (4.30) becomes 
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( ? m “ *0 ( f i-l ) + 0e,2 ( 7 i— 1 ) + <Pe,l (*i ) ~ <Pe, 2 (*i )]^ 7 7 1 


(4.31) 


Since the interpolation functions are expressed in terms of £ and T), they are combined with 
the transformed fundamental solution given by equations (4.18) and (4.19). The integrals of 
equations (4.29) are mapped into (g-ri) space to yield 

f" \tj ^^I dTd r.-i— r (,M) f'gfcn) (4.32a) 
Vl J dn j 0 dn 


f Jll^(r,r m - 0 ^^T«— jV(r(£),r7,r m -t 0 ) 

Jti-l J on ( ? m- ? 0) Jj 7(fi) J 0 


L r*rt f « Or 1 „,ir/ x df e (^,T]) dr] ( 4 . 32 b) 




The integrals are functions of £, 17 , the known boundary solution for t < r i _ 1 , the boundary 
conditions, and the unknown quantities at r = f,-. The expression does not admit analytic solution 
over space and time, but manipulation of the temporal integral is fruitful. 

The multiplication of the fundamental solution (eqs. (4.18) and (4.19)) and the interpolation 
function (eq. (4.31)) produces three core time integrals of the form 
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(4.33) 


The values of 77 at the limits of integration are found by applying its definition given by 
equation (4.17). Furthermore, the time integral is split into two parts to yield 
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(4.34) 


for k = 0, 1, 2 
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Both of the integrals on the right-hand side of equation (4.34) are of the same form, and a general 
solution appears to be possible. Yet, the time difference in the exponential function is (t m - 1 0 ) while 
the time difference in the numerator of the upper limit of integration is ( t m - /,•); this fact suggests 
that the integral must be computed for all combinations of ( t m - 7,) and (r m - fg) for 0 < i < m and 
1 < m < m max . Fortunately, the integrals can be transformed to yield standardized integration limits 
of 0 < 77 < 1. The results for each of the three types of integrals are 


JS 
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4 a(t m -t 0 )T] 


t m -tnJ 0 


A 


— r 


exp 
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(4.35a) 
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(4.35b) 
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(4.35c) 


Notice that the integration limits and the time difference in the exponential function have been 
standardized. Combined with the spatial integration, six component integrals result. The integrals are 
of standard form, and for a constant time step, are computed once, stored, and used repeatedly. The 
six integrals are denoted by the superscript tr to distinguish them from the steady-state integrals; 
they are given as 





(4.36a) 


(4.36b) 



(4.36c) 

(4.36d) 

(4.36e) 

(4.36f) 
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where 



4 a(t m -ti) 


(4.37) 


Furthermore, the time integration is performed analytically to yield 


'u=Jo 

(4.38a) 

I'{i = £|E 2 (^ 2 ©p| 

(4.38b) 


(4.38c) 

/&= 

(4.38d) 

,5 ’ i= W(l) eX ^ 2 ®^ 

(4.38e) 


(4.38f) 


where Ei and E 2 are exponential-integral functions whose expressions are given in reference 12. The 
integration over £ is performed numerically using Gaussian quadrature with the distance function 
given by equation (3.21). 

Summarizing up to this point, interpolation functions that are linear in both time and space are 
given for the dependent variables. The interpolation functions are expressed in terms of the trans- 
form coordinates £ and 77 . The multiplication of the interpolation functions to the kernels of the 
integral transform (i.e., the weight function and its derivative) yield three core temporal integrals as 
functions of r\ . The limits of integration and the time difference in the exponential term of the core 
integrals are standardized. These three core time integrals are combined with the spatial integration 
to produce six component integrals in space and time. The temporal integration is performed analyti- 
cally and the spatial integration is computed using Gaussian quadrature. Utilizing the definitions of 
the six integrals, the expression for the integral transforms of equations (4.29) are 
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Notice that the integrals Ijj appear in pairs of the form (q/^_ j -c 2 Ij r j) because of the splitting 
of the integration limits shown in equation (4.34). Also, terms of the form ( t m - tf) result from the 
standardization given by equation (4.35). The equations are tedious and the notation is abstruse, yet 
the concept is straightforward and should not be lost in the details of the algebra. The equations are 
obtained simply by multiplying the interpolation function by the transform kernels. 

Equations (4.39) are valid for any distribution of the time steps. Recall, though, that if values of 
t m are selected arbitrarily, then the integrals ifj must be calculated for all possible combinations of 
( t m - ti) for 0 < i < m and 1 <m< m max . If the time step is constant, however, then ( ) 
becomes At and (t m - tf) becomes an integer multiple of At. The integrals need to be calculated 
only times. Furthermore, because the influence coefficients are a function of a constant At, 

often only one matrix inversion is required for all time steps. Extracting out At, the integrals over 
time are 
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(4.40a) 
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4.4.2 Singular Integrals- The previous section develops expressions for the integral transforms 
when the integrand is nonsingular; this section presents expressions for the integral transforms for a 
singular integrand. The singularity occurs when the source point lies on the element defining the 
path of integration. Because the distance function simplifies in this orientation, the integral transform 
admits a complete analytic solution over time and space. Furthermore, since the source point lies 
along the path of integration ( r • n = 0), the integral of equation (4.29a) is zero. 

The singular integral requires an interpolation function for cTTg/dn. It is convenient to redefine 
the spatial transformation such that the origin lies on the source node point. The range of the trans- 
formation becomes -a < % < 1 - a , where a is the offset distance of the source point. The notation 
in this section is analogous to the notation used for Orientation 3 of the steady-state transform of 
Section 3.5.3. Letting subscript p denote the source point and q denote the other node point, the 
interpolation function for dt l e /dn is 



1 

(1 -b-a) 






dn 




(1 


dn 


(4.41) 


In the new transformation, the second node point is located at — 1 a b. 
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Similar to the derivation in the previous section, multiplication of the interpolation function to 
the transformed fundamental solution produces component integrals over time and space like the 
integrals of equations (4.36). With the revised spatial transformation, the component integrals are 
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(4.42a) 

(4.42b) 

(4.42c) 
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Because only one of the integral transforms is nonzero, only four component integrals are needed (as 
opposed to the six component integrals for the nonsingular case given by equations (4.36)). Each of 
the four integrals admits an analytic solution over time and space. The derivation of l[j is presented 
in detail as an example. The other three are derived similarly, but, for brevity, only the final expres- 
sions are given. With the distance function explicitly expressed in terms of t , , l{j is given by 
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The integral is singular at the origin. An additional transformation to polar coordinates allows 
the Cauchy principle value to be computed in the limit as the radius goes to zero. The domain is 
divided into two regions along the r\ axis. A separate polar coordinate transformation is prescribed 
for each region, as shown in figure 8. In each region, the limits of integration are easily determined 
by splitting the integral in two at the comers of %-rj space; i.e., at (£, 7?) = (1 - a , 1) and at 
(£ r?) = (-a, 1) . Both regions are transformed to polar coordinates and subdivided in two to 
yield four separate integrals 
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£=-a £ = 0 § = 1 -a 

Figure 8. Singular integration domain and polar transformation. 
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The inner integration over A is performed analytically to yield 
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(4.46) 


Each of the four integrals of the right-hand side of equation (4.46) is mapped to a dummy variable, 
u, through separate substitutions. Letting u = tan 0/(1 - a), u = tan 0 , u = cot y / , and u = cot y/ja for 
the first through fourth integral terms, respectively, yields 
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The first and fourth integrals on the right-hand side of equation (4.47) are of the same form, and the 
second and third integrals are of the same form. The solutions of the third and fourth integrals are 
addressed with the knowledge that the solutions of the first and second are obtained by substitution 
of (1 - a) for a. The fourth integral of equation (4.47) yields to analytic solution. Furthermore, the 
solution is expanded into series form; the result is 
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The third integral does not readily submit to analytic solution. Nevertheless, the integrand is 
expanded into series form and each term is integrable. The Cauchy principle value is found for 
each term at the upper limit of integration; the result is 
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The two series expansions found in equations (4.48) and (4.49) are combined to yield 
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Since the two series solutions can be combined, the sum of the solutions of the third and fourth 
integrals renders two terms: a series expansion given by equation (4.50) and an exponential integral 
function. The first and second integrals combine in like fashion. Summing all four integrals, the 
solution for l{j is 
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(4.51c) 
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The series shown in equations (4.51) are convergent for ft < 1. In practice, the element is broken 
into subelements such that the length of the subelement containing the singularity yields ft = 0.8 . 
For ft = 0.8, the series converges to sufficient accuracy in six terms. The remaining subelements are 
nonsingular and are integrated numerically as in the previous section. 

When a = 0, the expressions of equations (4.51) are undefined; however, the limiting form is 
found to be 
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The analytic expressions given by equations (4.51) and (4.52) are compared to numerical integration 
in figures 9 through 12. The numerical integration is performed using two-dimensional Simpson’s 
Rule with 10 9 functional evaluations. The analytic expressions are seen to reproduce the numerical 
values over a range of ft and a. Recall that the analytic expressions are not employed for ft > 0.8 
to ensure convergence of the series solution. The only discrepancy between the analytical and 
numerical solutions is for at small values of ft, as shown in figure 11. This difference is caused 
by error in the numerical solution; the Simpson integration requires extremely fine spacing to 
resolve the singularity of the integrand. 
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Figure 9. Comparison of analytical to numerical solution of transient, singular, component 
integral l[j. 



Figure 10. Comparison of analytical to numerical solution of transient, singular, component 
integral 
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Figure 11. Comparison of analytical to numerical solution of transient, singular, component 
integral I l {\. 



Figure 12. Comparison of analytical to numerical solution of transient, singular, component 
integral I'f , . 
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With I t ] r l defined, the expression for the singular integral transform is 
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4.5 Numerical Results 

The transient algorithm previously developed is tested against three analytic solutions. All three 
test cases are simple one- and two-dimensional geometries for which transient analytic solutions are 
found in reference 13. The first case is the temperature decay in a one-dimensional rod with an initial 
linear temperature distribution. The second case is the temperature decay of a two-dimensional 
plate initially at uniform temperature. Finally, the third case is a one-dimensional rod with a time- 
dependent fluctuating temperature boundary condition. 

4.5.1 One-Dimensional Rod- The first test case is a one-dimensional rod of length / with an initial 
linear temperature profile. The initial and boundary conditions are given by 

Initial conditions: T(x,t = 0 ) = (/ - x)/l for 0 < x < l 

Boundary conditions: dT/dn(x = 0 ,t) = 0 for t > 0 

T(x = 1,0 = 0 for t>0 

The domain is discretized using 1 1 elements along the length of the rod and 5 elements of length 0.1/ 
in the second dimension. Symmetry boundary conditions are imposed along y = 0 and y = 0.51 to 
simulate a one-dimensional problem. Temperature values at five select points along the length of the 
domain are compared to the analytic solution in figure 13. The abscissa is time nondimensionalized 
by the thermal diffusivity, a, and the length of the rod. The boundary element solution reproduces 
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Figure 13. Comparison of analytical to boundary element transient temperature profiles in one- 
dimensional rod. 


the analytic solution throughout the domain and for all time. The error, expressed as a percentage of 
the difference between the analytical and numerical solutions, is plotted in figure 14. The greatest 
error occurs at the initial time step; the cause of this effect is discussed in the next section. After the 
initial time step, the error decays to within 0.3 percent of the analytical solution. 

4.5.2 Two-Dimensional Plate- The second test case is a finite square with unit initial temperature. 
At time zero, the boundary temperature is instantaneously dropped to zero. Heat then flows out of 
the domain and the temperature drops to zero in the steady state. For this simulation, the boundary 
is discretized using 1 1 equally spaced elements per side. The square ranges from (x,y) = (0,0) to 
(jc ,y) = (1,1). Symmetry boundary conditions are imposed along x = 0 and y = 0 . The temperature 
is set to zero on the remaining two sides. In figure 15, the boundary element solution is compared 
to the analytic solution at five locations inside the domain. The points that are plotted, although 
selected for broad representation of the solution, are completely arbitrary and have no correspon- 
dence to any type of interior discretization. Any other five points could easily have been selected. 
Except for the first two or three time steps, the boundary element and analytic solutions are nearly 
identical to plotting accuracy. The error of the solution is presented in figure 16. Similar to the 
results presented in figure 14, the greatest error occurs at the initial time step. The error is slowly 
rising in time past 0.2 because the error is expressed as a percentage of the analytical solution, which 
is approaching zero. Thus, the percentage error is increasing even though the absolute error is 
decreasing. 
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at / 1 2 


Figure 14. Error of boundary element temperature profiles in a one-dimensional rod. 



at U 2 


Figure 15. Comparison of analytical to boundary element transient temperature profiles in a two- 
dimensional plate. 
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Figure 16. Error of boundary element temperature profiles in a two-dimensional plate. 


The slight deviations near time zero are explained by recognizing that the heat flux at time zero 
is infinite because of the instantaneous temperature drop. Consequently, the initial temperature 
gradient is computationally undefined which undermines the ability of the linear shape function to 
describe the thermal gradient over the first time step. This computational obstacle is circumvented 
by dropping the order of the temporal shape function over the first time step. Thus, the first time step 
for the transient calculation employs a constant temporal shape function while all subsequent time 
steps use a linear. As a result, the first time step has a reduction in accuracy as compared to the 
remainder of the computation. The effects of the accuracy reduction, however, do not propagate in 
time. 

4.5.3 One-Dimensional Rod with Fluctuating Boundary Condition- The final test case is a one- 
dimensional slab of length / with a time-dependent, fluctuating boundary condition. The initial and 
boundary conditions are given by. 

Initial conditions: T(x,t = 0) = 0 for 0 < jc < 1 

Boundary conditions: T{ x = 0,f) = 0 for t > 0 

T(x = l,t) = sin(10;rr + n / 2) for t > 0 
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The domain is discretized using 10 equally spaced elements on each side along the length of the 
domain plus 4 additional elements in the second dimension. The boundary element solution is 
compared to the analytic solution at 3 locations along the length of the domain in figure 17. 
Once again, the boundary element solution is identical to the analytic solution within plotting 
accuracy. 



Figure 17. Comparison of analytical to boundary element transient temperature profiles in one- 
dimensional rod with fluctuating temperature boundary conditions. 


The only marked error occurs at the first time step because of the constant temporal discretization 
mentioned previously. The L 2 error over all space is shown in figure 18 as a function of time and the 
time step. The error fluctuates at the same frequency as the boundary condition forcing function. 

This result is consistent with finite difference analysis, which shows the truncation error to posses 
the same growth property as the exact solution. (See ref. 14.) In this case, the temperature field is 
continually perturbed by the boundary condition, and, therefore, the solution and the error do not 
decay in time. Furthermore, the error decreases as the time step is decreased. The values of the peaks 
of the fluctuating errors are plotted as functions of the time step in figure 19. The slope of the curve 
shows the algorithm to be second-order accurate in time. Moreover, since the same interpolation 
function as the steady-state algorithm is employed, the algorithm retains second-order spatial 
accuracy. 
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Figure 18. Error of boundary element transient temperature profiles in one-dimensional rod with 
fluctuating temperature boundary conditions. 



0.001 0.01 0.1 1 

Time Step 


Figure 19. Error of boundary element transient temperature profiles in one-dimensional rod with 
fluctuating temperature boundary conditions. 
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5 Conclusion 


Two-dimensional, steady-state, and transient boundary element algorithms are designed for 
coupling with CFD flow solvers. The algorithms feature linear boundary elements for efficient 
coupling with discretized fluid domains that typically assume linear segments between grid points. 
Furthermore, the solution variables of land dT/dn are approximated using a linear interpolation 
function with offset node points; the offset nodes are implemented to uniquely define the nodal, 
element, surface normal. 

Analytic expressions for the requisite boundary integral transforms are derived for the steady- 
state algorithm. The analytic solutions are possible for both singular and nonsingular integrands 
because of the simplifications that result from using isoparametric linear shape functions. The 
steady-state algorithm is used to compute a variety of test cases that admit analytic solutions. The 
boundary element algorithm reproduces the analytic temperature to second-order spatial accuracy. 

The transient algorithm incorporates the transient fundamental solution approach. The boundary- 
only character of the solution procedure is maintained by originating the integral transforms from the 
initial time. The dependent variable interpolation functions are linear in both space and time. A 
series solution of the singular integral transform is derived and verified by comparison to numerical 
integration. The remaining integral transforms are computed using Gaussian quadrature. The 
transient algorithm is shown to be second-order accurate in the computation of test cases for which 
analytic solutions exist. 

The verification tests presented in this publication indicate that the BEM algorithms developed 
herein are suitable for accurate modeling of linear heat conduction. The algorithms are tailored for 
efficient and conservative coupling to CFD discretizations through the use of linear elements. 

Further modifications to the codes may be needed to incorporate boundary conditions required for 
coupling to a CFD code. 
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